Chapter 5 

Initial Value Problems 



5.1 Finite Difference Methods 

We don't plan to study highly complicated nonlinear differential equations. Our first 
goal is to see why a difference method is successful (or not). The crucial questions 
of stability and accuracy can be clearly understood for linear equations. Then we 
can construct difference approximations to a great variety of practical problems. 

Another property we need is computational speed. That depends partly on 
the complexity of the equation u' = f{u,t). Often we measure speed by the number 
of times that f{u,t) has to be computed in each time step (that number could be 
one). When we turn to implicit methods and predictor- corrector methods, to improve 
stability, the cost per step goes up but we gain speed with a larger step At. 

This chapter begins with basic methods (forward Euler, backward Euler) and then 
improves. Each time, we test stability on u' = au. When a is negative. At is often 
limited by —a At < C. This has an immediate effect: the equation with a = —99 
requires a much smaller At than a = — 1. Let me organize the equations as scalar 
and vector, nonlinear in general or linear with constant coefficients a and A: 

1 equation u' = f{u,t) u' = au a ~ df/du a < 

A?" equations ul = fi{u,t) u' = Au Aij'^dfi/duj ReA(A) < 

For an ordinary differential equation u' = f{u,t), good codes will increase the 
accuracy (and keep stability) far beyond the 0(At) error in Euler 's methods. You 
can rely on freely available software like ode45 to make the two crucial decisions: 

1. to choose an accurate difference method (and change the formula adaptively) 

2. to choose a stable time step (and change At adaptively). 

We will introduce accurate methods, and find the stability limits on At. 
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Stiff Differential Equations 



First we call attention to one extra question: Is the equation stiff ? Let me begin 
with a made-up example to introduce stiffness and its effects: 



e * + e 



T 



T 



controls decay controls At 



The step At is 99 times 
smaller because of e~^^* 
that disappears anyway 



Those decay rates —1 and —99 could be the eigenvalues of A, as in Example 1. 



Example 1 



d 
dt 
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The solution has v{t) = e * + e and w{t) = e' 



-99t 



The time scales are different 



by a factor of 99 (the condition number of A). The solution will decay at the slow 
time scale of e~*, but computing e~^^* may require a very small At for stability. 

It is frustrating to have At controlled by the component that is decaying so fast. 

Any explicit method will have a requirement 99At < C . We will see how this happens 
and how an implicit method (like odel5s and od23t in MATLAB) can avoid it. 

Trefethen [ ] points to these applications where stiffness comes with the problem: 

1. Chemical kinetics (reactions go at very different speeds) 

2. Control theory (probably the largest application area for MATLAB) 

3. Circuit simulation (components respond at widely different time scales) 

4. Method of Lines (large systems come from partial differential equations). 



Example 2 The N by N second difference matrix K produces a large stiff system: 

Method or Lmes — = has — — = — . 2 

dt (Ax)2 dt (Ax)2 ^ ' 

This comes from the heat equation du/dt = d'^u/dx^, by discretizing only the space 
derivative. Example 1 had eigenvalues —1 and —99, while Example 2 has N eigenvalues — 
but the difficulty is essentially the same I The most negative eigenvalue here is about 
a = — 4/(Ax)^. So a small Ax (for accuracy) will require a very small At (for stability). 

This "semidiscrete" method of lines is an important idea. Discretizing the 
space variables first produces a large system that can be given to an ODE solver. (We 
have ordinary differential equations in time.) If it wants to, that solver can vary the time 
step At and even the discretization formula as u(t) speeds up or slows down. 
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This method splits the approximation of a PDE into two parts. Finite differences/finite 
elements in earlier chapters produce the first part (discrete in space). The upcoming 
stability-accuracy analysis applies to the second part (discrete in time). This idea is very 
simple and useful, even if it misses the good methods later in this chapter that take full 
advantage of space-time. For the heat equation Ut = u^x, the useful fact Uu = u^xxx 
allows us to cancel space errors with time errors — which we won't notice when they are 
separated in the semidiscrete method of lines. 



Forward Euler and Backward Euler 

The equation u' = f{u,t) starts from an initial value u{0). The key point is that the 
rate of change u' is determined by the current state u at any moment t. This model 
of reality, where all the history is contained in the current state u{t), is a tremendous 
success throughout science and engineering. (It makes calculus almost as important 
as linear algebra.) But for a computer, continuous time has to change to discrete 
time. One differential equation allows many difference equations! 

The simplest method (Euler is pronounced "Oiler") uses a forward difference: 



Forward Euler 



Tl-l-l 



At 



fiUn, tn) is Un+1 = f/„ + At /„ 



(3) 



Over each At interval, the slope of U doesn't change. Figure 5.1 shows how the correct 
solution to u' = au follows a smooth curve, while U{t) is only piecewise linear. A 
better method (higher accuracy) will stay much closer to the curve by using more 
information than the one slope /„ = f{Un,tn) at the start of the step, 
replacements 




> t 




> t 



Figure 5.1: Forward Euler and Backward Euler for u' = u. One-step errors ~ |(Ai(:)^. 



Backward Euler comes from using fn+i at the end of the step, when t = tn+i^ 



Backward Euler 



n-l-l 



At 



= f{Un+l, tn+l) is Un+1 " At /„+i = f/„ . (4) 



This is an implicit method. To find f/„+i, we have to solve equation (4). When 
/ is linear in u, we are solving a linear system at each step (and the cost is low if 
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the matrix is tridiagonal, like I — At K in Example 2). We will comment later on 
iterations like Newton's method or predictor-corrector in the nonlinear case. 

The first example to study is the linear scalar equation u' = au. Compare forward 
and backward Euler, for one step and for n steps: 

Forward Un+i = (1 + a At)f/„ leads to U"^ = (1 + a At)"t/o • (5) 

Backward (1 - a At)Un+i = f/„ leads to = {1 - a At) '"^Uq . (6) 

Forward Euler is like compound interest, at the rate a. Each step starts with Un and 
adds the interest aAtUn- As the stepsize At goes to zero and we need T/ At steps to 
reach time T, this discrete compounding approaches continuous compounding. The 
discrete solution Un approaches the continuous solution of u' = au: 

{1 + aAt'f/'^^ approaches e"'^ as At ^ . 

This is the convergence oi U to u that we will prove below, more generally. It holds 
for backward Euler too, because (1 — a At)~^ = 1 + aAt + higher order terms. 

The stability question arises for very negative a. The true solution e~"*M(0) 
is extremely stable, approaching zero. (If this is your own money, you should change 
banks and get a > 0.) Backward Euler will be stable because it divides by 1 — a At 
(which is greater than 1 when a is negative). Forward Euler will explode if 1 + a At 
is smaller than —1, because its powers grow exponentially: 

Instability l + aAt<-l when aAt < -2 • i \ i (7) 

1+aAt -1 1 



That is a pretty sharp borderline between stability and instability, at —aAt = 2. For 
u' = —20m which has a = —20, the borderline is At = ^ = ^. Compare the results 
at time T = 2 from At = yt (22 steps) and At = | (18 steps): 

Stable = ^ (1 + a At)22 = (^-^^ ~ -012 

1 / 11\'' 

Unstable At = - (1 + a At)^^ = -— ^ 37.043 
9 V 9 / 

I would describe backward Euler as absolutely stable (A-stable) because it is stable 
whenever Re a < 0. Only an implicit method can be A-stable. Forward Euler is a 
stable method(!) because it succeeds as At — > 0. For small enough At, it is on the 
stable side of the borderline. 

In this example a good quality approximation requires more than stability (even 
At = is too big). Those powers of — alternate in sign, while e~^°* stays positive. 



Accuracy and Convergence 



Since forward and backward differences are first order accurate, it is no surprise that 
the errors from forward and backward Euler are 0(At). This error e = u — U is 
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replacements 




Figure 5.2: Forward Euler (stable and unstable, with ~aAt going above 2). 



measured at a fixed time T. As At decreases, so does the new error added at each 
step. But the number of steps to reach that time increases, keeping nAt = T. 

To see why the error u{T) — U{T) is 0{At), the key is stability. We need to know 
that earlier errors don't increase as they are carried forward to time T. Forward Euler 
is the simplest difference equation, so it is the perfect example to follow through first. 
(The next sections will apply the same idea to partial differential equations.) Euler 's 
difference equation for u' = f{u,t) = au is 

Un+l = Un + At f{Un, t„) = f/„ + « At f/„ . (8) 

The true solution at time nAt satisfies (8) except for a discretization error DE: 

Un+l = Un + AtUn + DE = Un + aAtUn + DE„+i . (9) 

That error DE is of order (At)^ because the second-order term is missing (it 
should be |(At)^M^, but Euler didn't keep it). Subtracting (8) from (9) gives a 
difference equation for the error e = u — U, propagating forward in time: 

Error equation e„+i = e„ + a Ate„ + DE„_|_i . (10) 

You could think of this one-step error DE.„ deposit like (At)^ into your account. 

Once the deposit is made, it will grow or decay according to the error equation. To 
reach time T = N At, each error DE^ at step k has N — k more steps to go: 

cn = {I + a Atf-^DEi + ■■■ + {l + aAt)^-''I)Ek + ■ ■ • + DEjv . (11) 

Now stability plays its part. If a is negative, those powers are all less than 1 (in 
absolute value) — provided 1 + aAt does not go below —1. If a is positive, those 
powers of 1 + a At are all less than (e"'^*)^ = e""^. Then the error ejy has terms 
in (11), and every term is less than c(At)^ for a fixed constant c: 

Error bound {snI = \un - Un\ < N c (Aty = cT At . (12) 

The errors along the way, of size (At)^, combined after = T/At steps into an 
overall At error. This depended on stability — the local errors didn't explode. 
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Notice that the error growth follows the difference equation in (10), not the dif- 
ferential equation. The stiff example with aAt= (— 20)(|) gave a large 1 + a At, 
even when e"^* was small. We still call forward Euler a stable method, because 
as soon as At is small enough the danger has passed. Backward Euler also gives 
I Cat I = 0{At). The problem with these first-order methods is their low accuracy. 

The local discretization error DE tells us the accuracy. For Euler that error 
DE ?a ^{At)^u" is the first term that Euler misses in the Taylor series for u(t + At). 
Better methods will capture that term exactly, and miss on a higher-order term. The 
discretization error (we find it by comparing u{t+At) with f/^+i, assuming u{t) agrees 
with Un) begins with a power (At)^"*""*^: 




The method decides c and p + 1. With stability, T/At steps give a global error of 
order [Aty. The derivative of u shows whether the problem is hard or easy. 

Error estimates like (13) appear everywhere in numerical analysis. The 1, —2, 1 
second difference has error y^{Ax)*u"" . Partial differential equations (Laplace, wave, 
heat) produce similar terms. In one line we find c = — | for backward Euler: 



[Un+l - Ur 



^tuu, - (At<+^^«:) -At«+Ato ^ (14) 



The global estimate \u — U\ < C At max|M"| shows no error when u is linear and 
u" is zero (Euler's approximation of constant slope becomes true). For nonlinear 
equations, the key is in the subtraction of (8) from (9). A "Lipschitz" bound L on 
df /du replaces the number a in the error equation: 

\fiu,t)- fiU,t)\<L\u-U\ gives e„+i < (1 + LAt) e„ + DE„+i . (15) 



Second-Order Methods 

Here is a first way to increase the accuracy. We could center the equation at the 
midpoint {n + |)At of the step, by averaging f{Un,tn) and /(f/n+i, Wi)- The result 
is an implicit second-order method use in MATLAB's ode23t: 

Trapezoidal rule/Crank-Nicolson — — —{fn+i + fn) • (16) 

So Un+i — I At fn+i = Un + \ At fn- FoT OUT model u' = f{u) = au, this is 

(1 - la At) Un+i = (1 + At) Un which gives f/„,+i = . (17) 

2 2 \1 — ja At / 

The true solution has Un+i = e°^*M„. Problem 1 will find DE ^ (At)^. The equation 
is stable. So = T/At steps produce |eAr| = \un — Un\ < cT(At)^. 
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How to improve forward Euler and still keep it explicit ? We don't want f/„+i on 
the right side of the equation, but we are happy with f/„_i (from the previous step !). 
Here is a combination that gives second-order accuracy in a two-step method: 

"Adams-Bashforth" ^" = ^ /(C/„, U) - ^ /(C/r^-l, ^r^-l) • (18) 



Remember that /(f/„_i, t„_i) is already computed in the previous step, going from 
n — 1 to n. So this explicit multistep method requires no extra work, and improves 
the accuracy. To see how (At)^ comes out correctly, write u' for /: 

<-i -\<-\ ^< - <^ = + ^ ""'^ ■ 

Multiplied by At in (18), that new term |(At)^M'^ is exactly what Euler missed. Each 
extra term in the difference equation can increase the accuracy by 1. 

A third possibility uses the already computed value t/„-i (instead of the slope 
fn-i)- With |, — |, I chosen for second-order accuracy, we get an implicit backward 
difference method: 

Backward differences 2^^t ~ f\Un+iitn+i) ■ (19) 



What about stability? The implicit method (17) is stable even in the stiff case, 
when a is very negative. 1 — |a At (left side) will be larger than 1 + |a At (right side). 
(19) is also stable. The explicit method (18) will be stable if At is small enough, but 
there is always a limit on At for stiff systems. 

Here is a quick way to find the stability limit in (18) when a is real. The limit 
occurs when the growth factor is exactly G = — 1. Set Un+i = — 1 and Un = 1 
and Un-i = —1 in (18). Solve for a when f(u,t) = au: 

—2 3 1 

Stability limit in (18) — = - a + - a gives a At = -1 . So C =1. (20) 



Those second-order methods (17)-(18)-(19) are definitely useful ! The reader might 
suggest including both t/^-i and to increase the accuracy to p = 3. Sadly, this 
method is violently unstable (Problem 4). We may use older values of U in backward 
differences or f{U) in Adams formulas, but including both U and f{U) for even higher 
accuracy produces instability for all At. 



Multistep Methods: Explicit and Implicit 

By using p older values of either U or f{U) (already computed at previous steps!), 
the accuracy can be increased to order p. Each VU is U{t) — U{t — At) and p = 2 
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is (19): 

Backward differences f V + ^ + ■ ■ ■ + - V^) f/„+i = At f(Un+i, t„+i) • (21) 

V 2 p / 

MATLAB's stiff code odel5s varies from p = 1 to p = 5 depending on the local error. 
Using older values of the right side f{U) gives an Adams-Bashforth method: 

Un+l-Un = At{bifn + --- + bpfn-p+l) with fn = f{Un,tn). (22) 

The table shows the numbers 6 up to p = 4, starting with Euler for p = 1. 



order of 
accuracy 


bi 


b2 




b. 


limit on aAt 
for stability 


constant c 
in error DE 


p = 1 


1 








-2 


1/2 


p = 2 


3/2 


-1/2 






-1 


5/12 


p = 3 


23/12 


-16/12 


5/12 




-6/11 


3/8 


p = A 


55/24 


-59/24 


37/24 


-9/24 


-3/10 


251/720 



The fourth-order method is often a good choice, although astronomers go above p = 8. 
At the stability limit G = -1 as in (20). The local error DE ^ c(At)P+in(P+^) is a 
problem of step control. Whether it is amplified by later steps is a problem of stability 
control. 

Implicit methods have an extra term co/n+i at the new time level. Properly 
chosen, that adds one extra order of accuracy — as it did for the trapezoidal rule, 
which is the second method in the new table. Backward Euler is the first: 



order of 
accuracy 


Co 


Cl 


C2 


C3 


limit on aAt 
for stability 


constant c 
in error DE 


p = 1 


1 








— oo 


-1/2 


p = 2 


1/2 


1/2 






— oo 


-1/12 


p = 3 


5/12 


8/12 


-1/12 




-6 


-1/24 


p = 4 


9/24 


19/24 


-5/24 


1/24 


-3 


-19/720 



Every row of both tables adds to 1, so m' = 1 is solved exactly. 

You see that the error constants and stability are all in favor of implicit methods. 
So is the user, except when solving for f/„+i becomes expensive. Then there is a 
simple and successful predictor- corrector method, or else Newton's method. 



Predictor- 
Corrector 



Use the explicit formula to predict a new f/*+i 
Use to evaluate the right side /* 



''n+l 



n+l 



Use /„ , ]^ in the implicit formula to correct to a new Un+i- 
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The stability is much improved if another E step evaluates fn+i with the corrected 
Un+i- In principle we should continue the corrector to convergence. Frequently 1-2 
corrections are enough, and the extra calculation can be put to use. By comparing the 
predicted f/^+i and corrected f/„+i the code can estimate the local error and change 
At: 

local error DE ^ ^^^{Un+i - f/*+i), (23) 

where c* and c are the error constants in the tables for the predictor and corrector. 

Implicit methods often have a Newton loop inside each time step, to compute 
Un+i- The kth iteration in the Newton loop is a linear system, with the Jacobian 
matrix A\j^ = dfi/duj evaluated at the latest approximation U^_^^ with t = 

Here is the kth i to solve Un+i — Co/(?7„+i, = old values: 

Newton iteration (/ - coA('=))(t/i^V^ - = ('of{ui%, t„+i) + old values .(24) 

When f{u) = Au is linear, you only need one iteration (fc = starts with U^}_i = Un). 
For nonlinear f{u), Newton's rapid convergence squares the error at each step (when 
it gets close). The price is a new evaluation of f{u) and its matrix A^''^ of derivatives. 



Runge-Kutta Methods 

A-stable methods have no stability limit. Multistep methods achieve high accuracy 
with one or two evaluations of /. If more evaluations are not too expensive, Runge- 
Kutta methods are definitely competitive (and these are self-starting). They are 
compound one-step methods, using Euler's Un + At fn inside f : 

Simplified Runge-Kutta ^" = ^ [/„ + /(f/„ + At /„, t„+i)] . (25) 

You see the compounding of /. For the model equation u' = au the result is 

Un+l = Un+ -At\aUn + Ci{Un + At aUn)\ = (l + CtAt + -O^ At^^Un = G Un ■ (26) 

This confirms the second-order accuracy; the growth factor G agrees with e"^* through 
(At)^. There will be a stability threshold a At = C, where |G| = 1: 

Stability limit IG] = 1 1 + C + ^ | = 1 (for complex C = a^ih) 

The complete stability limit is not a point but a closed curve in the complex plane. 
Figure 5. shows all the complex numbers C at which |G| = 1. 
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The famous version of Runge-Kutta is compounded four times and achieves p = A: 

Runge-Kutta Un+i - Un = + + + k^) (27) 

/^t 3 



^1 = 2 ^ 2 

For this one-step method, no special starting instructions are necessary. It is simple 
to change At as you go. The growth factor reproduces e"^^* through ^ a'^At^. The 
error constant is the next coefficient Among highly accurate methods, Runge- 
Kutta is especially easy to code and run — probably the easiest there is. MATLAB's 
workhorse is ode45. 

To prove that the stability threshold a At = —2.78 is genuine, we reproduce the 
solution of m' = —100m + 100 sin t. With u{0) = 0, Runge-Kutta gives 

Ui2Q = 0.151 = m(3) with At 
f/ioo = 670,000,000,000 with At 



120 
3 

Too 



and aAt = —2.5 
and aAt = — 3 



(2^ 



Differential-Algebraic Equations 

Our system of equations might not come in the form u' = f{u,t). The more general 
form F{u',u,t) = appears in many applications. It may not be easy (or possible) 
to solve this system explicitly for u'. For large electrical networks (VLSI on chips), 
transistors produce a highly nonlinear dependence on u' and u. 

Here is a simple model in which solving for the vector u' is impossible: 

Mu' = f{u,t) with a singular matrix M (rank r < A^) . (29) 

Suppose we change variables from u to v, so as to diagonalize M. In the new variables 
we have r true differential equations and N ~ r algebraic equations {no derivatives) . 
The system becomes a differential- algebraic equation (DAE): 

dvi +\ ■ 1 

DAE dt 1, , TV, ; ^^^^ 

= fi{vi,. . . ,VN,t) i = r + l,...,N. 

The N — r algebraic equations restrict f to a surface in A^-dimensional space. The r 
differential equations move the vector v{t) along that surface. 

Petzold [46] developed the DASSL package for solving DAE's in the general form 
F[u', u, t) = 0, replacing u' by a backward difference of u. The best recommendation 
we can make is to experiment with this software ! 
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Problem Set 5.1 

1 The error in the trapezoidal (Crank-Nicolson) method ( ) comes from the 
difference 

/ aAt\ / a At /aAt\2 \ 

+ — + — + (— ) +■■■). 

This involves the two series that everyone should learn: the exponential series 
for e"^* and the geometric series 1 + x + + ■ ■ ■ for . 

^ l—x 

Multiply inside the brackets to produce the correct |(aAt)^. Show that the 
(aAt)^ term is wrong by c = ^. Then the error is DE ^ j^{AtY u'" . 

2 Try Runge-Kutta on u' = -lOOn + lOOsint with At = -.0275 and -.028. 
Those are close to the stability limit —.0278. 

3 For the backward difference error in (19), expand |(3e"'^* — 4 + e~''^*) — aAt e"'^* 
into a series in aAt. Show that the leading error is —^{aAt)^ so that c = — |. 

4 Stability in (19). 



^aAt 



1 + (aAt/2) 
1 - (aAt/2) 



